function y = func_ns(q, zeta, m, CONSTS)

    k0 = CONSTS.k0;
    a  = CONSTS.a;
    c  = CONSTS.c;

    y_vec =  func_total(q, zeta, m, CONSTS); %func_total_new(q, zeta, m, CONSTS); 
    y_s_vec = - ((2i*pi*k0^2*a)/c) * func_1(q, zeta, m, CONSTS) + ...
                       (2i*pi/(c*a)) * func_2_for_ns(q, zeta, m, CONSTS);
    y = y_vec - (y_s_vec); 
    
    if(false)
        figure; plot(q, imag(y_vec), 'b-', q, imag(y_s_vec), 'r-', q, imag(y), 'k');
%         figure; plot(q, real(y_vec), 'b-', q, real(y_s_vec), 'r-', q, real(y), 'k');
    end

end

function y = func_2_for_ns(q, zeta, m, CONSTS)

    k0    = CONSTS.k0;
    a     = CONSTS.a;
    eps   = CONSTS.eps;
    eta = CONSTS.eta;
    sigma = CONSTS.sigma;
    
    if isempty(q)
        y = [];
        return;
    end
    
    y = (m^2) * ...
        besselmx(double('J'), m, k0*a*(1/sigma)*q, 0).^2 .* ...
        func_V(((1/sigma)*q), m, CONSTS) .* ...
        (1/(sqrt(eps*eta))).*(sqrt(eta/eps)) .* ...
        exp(-k0*q*abs(zeta));
    
end

